Particle-Number Projection and the Density Functional Theory 
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I. INTRODUCTION 

Superconductivity plays a central role in describing 
low-temperature properties of correlated many-fermion 
systems. Within the mean-field theory, fermionic pairing 
is best treated in the Hartree-Fock-Bogoliubov (HFB) 
[l| or Bogoliubov-de Gennes (BdG) [H formalism. In 
the presence of superconducting condensate, the stan- 
dard product state ansatz for the nuclear wave function 
breaks the particle- number (PN) symmetry [H, In 
principle, the broken symmetry needs to be restored, es- 
pecially if one looks at observables that strongly depend 
on PN. The many-body correlations associated with the 
symmetry-breaking are particularly important for small 
systems where the finite-size effects are appreciable, such 
as atomic nuclei or metallic grains, or in the limit of weak 
pairing where pairing correlations have dynamic charac- 
ter. 

For complex superconducting system, s a theoretical 
tool of choice is the Density Functional Theory (DFT) 
|, i]. The theory is built on theorems showing the 
existence of energy functionals for many-body systems, 
which include, in principle, all many-body correlations. 
The generalization of the DFT to the case of fermionic 
pairing was formulated for electronic superconductors in 
Refs. [sSH]- The resulting HFB/BdG equations can be 
viewed as the generalized Kohn-Sham equations of the 
standard DFT. 

In the nuclear case, the DFT is the only tractable 
theory that can be applied across the entire table of 
nuclides. Historically, the first nuclear energy density 
functionals appeared long ago [H, [l(| EH in the context 
of the Hartree-Fock (HF) method used with zero-range, 
density-dependent interactions such as the Skyrme force. 
The main ingredient of the nuclear DFT [12| is the en- 
ergy density functional that depends on densities and 



currents representing distributions of nucleonic matter, 
spins, momentum, and kinetic energy, as well as their 
derivatives (gradient terms). To account for nuclear su- 
perfluidity, the functional is augmented by the pairing 
term (see Ref. [l3( for a review). The challenges faced 
by the nuclear DFT are: (i) the existence of two kinds 
of fermions; (ii) the essential role of pairing; and (iii) the 
need for symmetry restoration in finite, self-bound sys- 
tems. The two latter points are of particular importance 
in the context of this study. The features (i) and (iii) are 
specifically nuclear; with very few exceptions, they are 
not present in the electronic Coulomb problem. 

It is important to recall that the realistic energy den- 
sity functional does not have to be related to any given 
effective Hamiltonian, i.e., an effective interaction could 
be secondary to the functional. This strategy is used 
in all modern nuclear DFT applications. In the absence 
of a Hamiltonian (and wave function), the restoration of 
spontaneously broken s ymm etries in DFT poses a con- 
ceptional dilemma [3, Il5l [T5 , [TtJ and a serious chal- 
lenge that needs to be properly addressed. One impor- 
tant question related to DFT for self-bound systems con- 
cerns the functional itself: how do you construct it in 
terms of intrinsic (body-fixed) densities? While it is pos- 
sible to formulate the Kohn-Sham procedure in language 
of intrinsic densities [H, dl, the pathway to practical 
applications is still not clear. 

Sticking to DFT for superconductors and PN symme- 
try, several schemes can be adopted. One is to formulate 
the theory in language of the usual (particle) density only, 
without explicitly invoking the anomalous (pair) de nsity 
that is at the heart of the PN symmetry violation [20L 1 2 1|] . 
Another strategy is to incorporate the PN restoration 
procedure into the DFT framework. This can be done 
by employing the generalized Wick's theorem (see, e.g., 
[22|, [23, [2J| ) ■ Recently, full PN projection before varia- 
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tion has been carried out for the first time within the 
Skyrme-DFT framework employing zero-range pairing 
|23l . \2SL l26l ] . It was demonstrated that the resulting pro- 
jected DFT equations (similar to the PN-conserving HFB 
equations originally proposed in Refs. [27|,|28|) can be ob- 
tained from the standard Skyrme-HFB equations in coor- 
dinate space by replacing the intrinsic densities and cur- 
rents by their gauge-angle-dependent counterparts. Us- 
ing the variation-after-projection method, one can prop- 
erly describe transitions between normal and supercon- 
ducting phases in finite systems, which are inherent in 
atomic nuclei. 

As mentioned above, the restoration of broken sym- 
metries in the framework of DFT causes a number of 
questions, mainly related to the density dependence of 
the underlying interaction and to different treatment of 
particle- hole and particle-particle channels [H, lipl. l29l . 
For instance, it has been realized for some time [22|, [25|, 
l29l l30l I3H |32| that the PN projection applied within 
the DFT framework is plagued with difficulties related to 
vanishing overlaps between gauge-rotated intrinsic states. 
This concerns any functional that uses density-dependent 
terms and thus is not related to an average of a Hamilto- 
nian. In particular, the most frequently used approaches 
based on the Skyrme, Gogny, or relativistic-mean-ficld 
functionals all fall into this category. 

In this study we investigate the analytic structure of 
the projected DFT, focusing on origins of difficulties. In 
recent works [33L l34j , a way to remedy some of the prob- 
lems has been proposed. The PN-projected Skyrme-DFT 
formalism employed in our work has been outlined in 
Ref. (23[, and we follow their notation. Our manuscript 
is organized as follows. The analytic structure of the pro- 
jected HFB is discussed in Sec. [TTJ The DFT extension 
of the formalism is described in lllll Numerical examples 
are contained in IIV1 Finally, Sec. fVl contains conclusions 
of this work. 



II. PARTICLE-NUMBER-PROJECTED HFB 

In the context of HFB theory [l| , the particle- number- 
projected (PNP) state is given by the standard expression 
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where N is the PN operator, 4> is the gauge angle, P/v is 
the projection operator for N particles, and |<I>) is the 
HFB wave function (generalized product state) which 
does not have well-defined particle number. This expres- 
sion, after the integral is discretized, is most often used 
in practical calculations. However, it only constitutes a 
specific realization of a more general form [35| given by 
the contour integral, 
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where C is an arbitrary closed contour encircling the ori- 
gin z = of the complex plane. 



A. Shifted HFB states 

Let us introduce several useful notations that will be 
used later. First, we call the operator appearing under 
the integral ([2]) the shift operator, 



- ,N _ Jv+i<P)N 



z(z) = z = e 



(3) 



parametrized by means of a single complex number z, 
ln(z)=?7 + The shift operator z(z) is parametrized 
by the complex number z and constitutes a non-unitary 
Bogoliubov transformation (in fact, a non-unitary single- 
particle basis transformation) of simple kind, i.e., 



or 



za^z * = za^, 



(4) 



(5) 



Obviously, for z = 1, the shift operator is equal to iden- 
tity. 

Second, we define the shifted HFB states as 



!*(*)) = *(*)|$>. 



(6) 



When the HFB state |$) is expressed through the Thou- 
lcss theorem Q (we assume an even number of particles 
for simplicity), 



|$)=AAexp(i^Z i ;„a+a+) |0), 

\ mn / 



(7) 



|0>, (8) 



the shifted HFB states read 

|$(z)) = N exp f \z 2 Z m n aUn 

\ mn J 

where M is the normalization constant of the HFB state 
([7]) . Similarly, for the HFB state expressed in the canon- 
ical basis or for a BCS state, 



I*) = II (u n + v n a+a+) |0), 



(9) 



n>0 
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the shifted state reads 

!*(*)) = II K + z 2 «»<44) |o>, (10) 

n>0 

where u n and v n are the real HFB occupation amplitudes 
in the canonical basis and the product Jln>o involves 
only one state from each pair of canonical partners (see 
Ref. Q for details). 

We call ziz) a shift, because it moves the HFB state 
|$) = |$(f)) from its original position at z = 1 to a dif- 
ferent point z in the complex plane. Since consecutive 
shift transformations correspond to products of the shift 
parameters z, the parameters r\ and <j> in Eq. ([3]) are ad- 
ditive. 



3 



B. Projected HFB states 



C. HFB sum rules 



The Thouless theorem (J?} allows us to express the HFB 
state |<E>) and shifted HFB state |3>(z)) as sums of com- 
ponents having different particle numbers, 



|$(z))=AAf> 2fe ^|0). 



fe=0 



(11) 



(12) 



, mn Z mn a m a t is the Thouless pair- 



where Z + = | J2 r , 
creation operator. It then trivially follows that the 
shift transformation does not change any of the particle- 
number-projected states (J5J, 



(Z+) N / 2 



(13) 



but only scales the coefficients in the sum of Eq. (fl"2"l) . 

Since the shifted states (| 12[) are manifestly analytical 
in z, all closed contours C in Eq. ([2]) give, by the Cauchy 
theorem, the same result. Among them, the integral in 
Eq. {1} simply corresponds to the unit circle. 

The analyticity of |$(.z)) results in a simple and elegant 
representation of the projected state: 



\9 N ) =Pn\&) = Resz 

z = 



-JV-1 



(14) 



Indeed, in the sum of Eq. p2|) , only the term with N = 
2k particles is multiplied by l/z and thus contributes to 
the residue at z = 0. This observation allowed Dietrich, 
Mang, and Pradal [36[ to formulate the so-called method 
of residues for calculating all kinds of matrix elements 
involving the projected state I* at). For example, the 
average HFB energy of the projected state can be written 
as a ratio of two residues: 



^HFB — 
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Res z-k- 1 

z = 
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The invariance of the projected state with respect to 
the integration contour can be formulated in another 
way; namely, one can utilize the property that an ar- 
bitrarily shifted HFB state can be equally well used to 
project the particle number. Indeed, for 



1 
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dzz N - N - 1 



we trivially have 



\* N (zo))=zS\H! N ), 



(16) 
(17) 



i.e., projection from a shifted HFB state changes only the 
phase and normalization of the projected state. We refer 
to this property as shift invariance. 



Since the HFB state (fTTj) is a superposition of projected 
states (TT3"1) . 



i$) = \*n) 



(18) 



N=0 



the HFB energy £hfb, 

£hfb = ($|#|$>, (19) 
can be expressed as the sum of projected energies (fl"5j) . 

oo 

•Ehfb = n) Ehfbi 

N=0 



(20) 



weighted by probabilities (*Ar|*Ar) = ($|*Ar) = ($|PAr|$) 
of finding a given PN component in the HFB state. 
Expression ([2H)l constitutes a useful sum-rule condi- 
tion, which has to be obeyed by any Hamiltonian-based 
HFB+PNP approach, and can be used to test the nu- 
merical precision of PNP techniques. 

A similar sum rule holds for any shifted state 



Mz )) = £ \* N (z )), 



N=0 



($(z )|ff|$(z )) = J2 |*>| aW <<M**>3 



N 

HFBi 



(21) 



(22) 



where the average energy of the shifted and unnormalized 
HFB state is related to its HFB energy -Ehfb^o) &s 



-^hfb(^o) 



(j>(z )|g|$(z )) 



(23) 



Finally, the sum rule for the non-diagonal matrix ele- 
ments can be written as: 



<$|H|$(z )) = J2 z^^n^n)^ 



N 

HFB- 



(24) 



N=0 



D. Transition matrix elements and transition 
densities 

Calculation of the matrix elements in Eq. (fT5")) between 
the original and shifted HFB states is straightforward, 
because the shifted states also belong to the family of the 
HFB states. In particular, their overlap is given by the 
Onishi formula [1], which in the canonical basis reduces 
to a simple expression, 



mm) = WK + z 2 vi). 



(25) 



n>0 
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Similarly, the generalized Wick's theorem [l[ can be used 
for evaluation of Hamiltonian matrix elements, 



($|ff|$(z)> = {$Mz))Ehfb(j>„Xz,Xz 



the pairing term impair that depends on the pairing tran- 
sition densities. It was first realized in Ref. [301 ], and 
then discussed by several authors [U [H, H3] , that the 

(26) 

poles are not cancelled separately in £fi e id and E pa i T , but 



where the so-called HFB transition energy density 
E-rfb(Pz,Xz,Xz) is a function of the shifted particle and 
pairing transition density matrices, 

Pz (m,rV) = (<S>\a+ a ,a ra \<P(z))/(<£Mz)) 

2 2 

X,(rtr,rV) = ($|(v CT ,a rCT |$(z))/($|$(z)) 



E 



Xz{vo,t'o') = ($|a+a+ CT ,|$(z))/($|$( 2 )) 



E 



rtT r'a' 

u n v n * r \n I it i\ 



The transition density matrices become the standard 
density matrices in the limit of z — ► 1. For simplicity, 
we do not explicitly show the isospin variables; this is 
not essential in the context of the present work. (See 
Ref. [HI for a complete formulation.) 



E. Poles of transition densities 

It is seen immediately from Eq. (f2"T|) that the transition 
density matrices have imaginary axis poles at 



±i\u n /v n \, 



(28) 



and, therefore, are not analytical. These poles carry over 
to the HFB transition energy density as well. The poles 
appear beyond the origin, z n ^0, provided all amplitudes 
u n are nonzero; we assume this hereafter, i.e., none of 
the canonical states is being blocked. We can also safely 
assume that all amplitudes v n are nonzero, because oth- 
erwise the corresponding states would not contribute to 
the density matrices at all. Of course, if there exist poles 
in the HFB transition energy density, they must be can- 
celled by the norm overlap ($|<f>(z)), because the Hamil- 
tonian matrix element (3>\H\Q(z)) is an analytical func- 
tion of z. 

However, as we discuss in the next section, whenever 
the transition energy density is not related to a Hamil- 
tonian, or some approximations are involved in Hamil- 
tonian's construction, the presence of the poles (f2"8")) re- 
quires special attention. For example, the exact HFB 
transition energy density, 

£hfb (p z , Xz , Xz ) = E kin (t z )+ ^fioid (p z ) + £ pair (xz , X.z ) , 

(29) 

is often split into the kinetic term -E'kin(Tz) that depends 
on the kinetic transition density, the mean-field term 
infield that depends on the particle transition density, and 



pair: 

only in the sum thereof, i.e., for the total HFB energy 
calculated for a given Hamiltonian. 

As the origin of the pairing interaction is believed to 
be different from that of the effective interaction in the 
particle-hole direction, it is customary to employ differ- 
ent Hamiltonians to calculate -Efidd and i? pa ir- This, how- 
ever, leads to a non-analytical behavior of -Ehfb due to 
the presence of poles in the complex z-plane; hence, to 
a priori contour-dependent projected HFB energies. We 
discuss this question in the next section in the more gen- 
eral context of the DFT energy functional. 



III. PARTICLE-NUMBER-PROJECTED DFT 

According to the DFT, the energy density of the sys- 
tem, £'dft(/9, X, X*)y can be written as a function of the 
local particle p(r) and pairing x( r ) densities obtained as 
the diagonal elements of the corresponding density ma- 
trices: 

X(r) = T,J- 2cr )x(ro-,r,~a) = J2 niJ u n v n \ip n (ra)\ 2 . 

(30) 

The nuclear density functionals for time-even systems 
also depend on kinetic r and spin-orbit J densities. An 
even larger set of densities enters the energy density for 
time-odd systems [l3l l38j . For simplicity, we discuss here 
the dependence on the particle density only, because ex- 
tension to other densities is straightforward. 

We note in passing that the densities corresponding to 
the shifted HFB state ^(z)) @ can be written as: 



p z (r) 
X 2 (r) 



E 

■n 

E 



l*IV 



< + M 4 ^ 



K + M K 



Ei^«( r<T )i 2 ' 

</?„(ra)| 2 . 



(31) 
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A. Transition energy density 

In the DFT approach, the Hamiltonian of the system 
does not appear explicitly; hence, the projected energy 
cannot be calculated as its expectation value in the pro- 
jected state. However, since the DFT energy density is 
most often postulated, not derived, we can apply the 
same philosophy to the projected energy, i.e., we can pos- 
tulate the projected functional. In doing so, we have to 
guarantee that it reverts to the projected HFB energy 
(p~5|) when the system is described by a Hamiltonian. In 
the present study, we do not discuss the construction of 
the projected DFT functional, but simply assume, as in 
most calculations up to now, that the DFT transition 
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energy density Euft(Pz,Xz,X.z) is the same as the DFT 
energy density -Edft(p, X, X*)j but with densities p, Xi 
and x* (ESS) replaced by the transition densities p z , Xz, 
and Xz (|2"T)) . This guarantees that in the limit of z — > 1, 
the projected functional gets back to the usual form. 

Since the overlap (|25|) and HFB transition energy den- 
sity (|26p depend only on the shift parameter z and not on 
its complex conjugation z*, it is natural to restrict fur- 
ther considerations to the DFT transition energy density 
parametrized in the same way, i.e., 



£'dft( z ) 



E, 



DFT 



CO- 



(32) 



Moreover, by construction, the DFT transition energy 
density depends only on z 2 , and therefore it must be a 
symmetric function of z, 



Edft(-z) — Euft(z). 



B. Projected DFT energy 



(33) 



Based on the above discussion, we postulate the pro- 
jected DFT energy in the form: 



N 

DFT 



§ c dzz~ N - 1 (<f | $(z))£ DFT (Pz, Xz , Xz ) 



2ni Res z 

z = 



-N-l, 



(34) 



At variance with the Hamiltonian-based HFB theory, 
the projected DFT energy may depend on the integra- 
tion contour C. Moreover, the numerator in Eq. (|3"4"1) 
is, in general, not equal to the residue at z—0 as in 
Eq. (fTS]) . Consequently, both the transition energy den- 
sity and the contour C define the projected energy in 
DFT. Since the projected DFT energy (13"4")) must be real, 
in view of condition (|32p , we restrict further consider- 
ations only to contours which are symmetric with re- 
spect to the real z-axis. Accordingly, only the upper-half 
contour C+ above the real axis can be considered and 
§ c dz--- = 2Re[§ c dz---]. 



C. Analytic properties 

In order to proceed further, we must investigate the 
analytic structure of the integrand £n(z) appearing in 
the numerator of Eq. 



£ N {z) 



„-N-l 



<£\<I>(z)}E UFT (p ZJ Xz,Xz 



(35) 



Let us first discuss the case when the DFT energy density 
Et>ft(Pi Xi X*) is a polynomial in local densities; hence, 
the DFT transition energy density £t>ft(Pz, Xzi Xz) is a 
polynomial in transition densities. The case of fractional 
power dependence requires special attention and will be 
discussed in Sec. IIII Fl 

Within the polynomial assumption, poles of the tran- 
sition densities (|28p do or do not appear as poles of the 
integrand (|3"5"|) . depending on the structure of the DFT 




FIG. 1: (color online) Schematic illustration of the analytic 
structure of the integrand (|35[) in the complex z-plane (see 
text). Small crossed circles denote imaginary axis poles. 
Three integration contours (CO, CI, C2), symmetric with re- 
spect to the real axis, are indicated. The poles having particle 
character (corresponding to canonical states lying above the 
chemical potential A) are located outside the unit circle CO 
while the hole poles lie inside CO. The unprojected ground 
state wave function corresponds to z=l. 



transition energy density. [For instance, quadratic (p=2) 
and cubic (p=3) terms are characteristic of two-body and 
three-body interactions, respectively] On the one hand, 
each polynomial term of the order p in densities (|30p 
brings about a pole of the order p. On the other hand, 
each term in the overlap (|28p produces a zero of the or- 
der q, where q is the degeneracy factor of the HFB den- 
sity matrix with the two-fold Kramers degeneracy not 
counted. (Note that the product in Eq. (f28|) contains 
only one term for each canonical pair.) In particular, for 
the spherical shell of angular momentum j, the degener- 
acy is q = j + \. 

When the poles of transition densities and zeros of 
the overlap ($|$(z)) are combined, the poles in £n{z) 
are of the order p — q. For single-particle states that 
have only two- fold Kramers degeneracy (q=l), and for 
the terms with p=2, one obtains the first-order poles in 
£n{z) with, in general, non-zero residues. Non- vanishing 
residues may also appear for higher-order poles corre- 
sponding to terms with p > 2. On the other hand, for 
four-fold degenerate states with q=2, terms with p=2 do 
not produce poles in £n(z), and only terms with p > 2 
may give rise to poles with non- vanishing residues. As 
discussed in detail in Ref. 30], for the energy density 
derived from a Hamiltonian, additional cancellations be- 
tween terms originating from particle-hole and pairing 
channels occur, and the first-order poles disappear. 

In Fig. [1] we schematically illustrate the analytic struc- 



6 



6 

Cl 



i 


i Im[z] 


^ 


y MMMB, 


















■■l 




1111111 














^\ 













Im[z] 







/ ^ 

H 


5 \ 
»4— — 








V < 


s> / 





Re[z] 



FIG. 2: Schematic illustration of analytic structure of the 
particle transition density at a fixed point in space r. The 
poles (crossed circles) and zeros (dots) of p z are located on 
the imaginary-.? axis. The regions of real negative p z are 
shaded. Three circular integration contours CO, Cl', and C2' 
are indicated. See text for details. 



ture of the integrand (J35J ■ Crossed circles on the imagi- 
nary axis represent poles of £n(z). Apart from the pole 
at z = 0, the integrand may have poles (|28[) distributed 
symmetrically in pairs with respect to the real axis. Poles 
located within the unit circle (CO in Fig. []} correspond 
to the canonical states with occupation numbers larger 
than 0.5, or with u n /v n < 1, i.e., with canonical energies 
below the Fermi energy A. Similarly, poles outside the 
unit circle correspond to canonical states lying above the 
Fermi energy. 

The unprojected HFB ground state |$), located at 
z=l, is shifted along the integration contour C, and its 
overlap and DFT transition energy contribute to the in- 
tegrand of Eq. ([55*)) . Standard projection formula ((T|) 
corresponds to the unit circle CO. Contours Cl and C2 
encircle a fewer number of poles in £jy (z) , with contour 
C2 surrounding only the single pole at the origin. Shapes 
of these contours are irrelevant, and only the points at 
which they cross the imaginary axis matter. For exam- 
ple, contours Cl and C2, shown in Fig. [TJ are equivalent 
to circular contours Cl' and C2' of Fig. [21 the latter 
being more practical in calculations. If the residues of 
the poles inside the unit circle are non-zero, the three 
integration contours shown in Fig. [T] may give different 
projected energies. Of course, contours including poles 
located outside the unit circle (not shown in Fig. [1]) may 
still give different results. 



D. Residues 

Let us now discuss the residues of the integrand (|35[) . 
From Eqs. (["2"5")) and ([3"3")1 . we see that the integrand is an 
odd function of z, 



£n(-z) = -£ N (z). 



(36) 



This is obvious for even particle numbers N, for which 
Eq. (fJU) has been derived, while for odd N, an additional 
power of z appears when shifting the blocked HFB state, 



|(j) odd ) 



no II ( u n + v n a+a^)\0), (37) 

n^in >0 



K d (z)> = za+ o J] (u n + z 2 v n a+a+)\0), (38) 



which gives 



($ odd |$ odd (z))=z -Q (u 2 

Tt^TtO>0 



~ 2 vl) 



(39) 



and renders the integrand ([33)) an odd function of z also 
for odd systems. 

Near the pole (|28j) . the term in the integrand that pro- 
duces the residue has the structure: 



£ N {z) 



n n {z) 

ui + z 2 v"i 



(40) 



where 1Z n (z) is an odd function of z, regular at the pole. 
Similarly, for the pole at z=0 we have 



£n(z) 



Ko(z) 



(41) 



Therefore, for pairs of poles that are symmetric with re- 
spect to z=0, the residues, 



(z^fi\Un I V n DTZn (z) 



Res £ N {z) = lim^ ±i i„ 

z — ±i\u n /v n I ■- 

Hn{i\u n /V n \) 

2i\u n v n | ' 

(42) 

have identical values. Hence, poles below and above the 
real axis yield the same contribution to the contour in- 
tegral. Based on this consideration, the projected DFT 
energy (|34p. expressed in terms of residues, reads: 



N 

DFT 



Res ($|$(.z)) 

2 = 



(43) 



(44) 



n=0 



where E^ FT (n) denotes the contribution from the nth 
pole, including the n=0 pole at the origin up to n = n 
(last pole encircled by C). 
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As an example, we explicitly calculate the residues for 
a term that depends on the squared particle density, 

Ev F t(Pz) = C J d 3 rp*(r), (45) 

with 

^w=E^ £ ^2 Ei^mi 2 (46) 



One can see that residues can be very large for poles 
corresponding to canonical states that have occupation 
numbers close to 1. These very large contributions to 
the projected DFT energy must be compensated by a 
similarly large contribution from the single pole at z=0. 
Therefore, within the DFT formalism, one cannot use 
the HFB expression (| 1 5 1) that involves only one residue 
at z=Q. 

Recall from our discussion in Sec. IIII CI that the poles 
have the order of p — q. In the above example, the poly- 
nomial order is p=2; hence, the residue (|T7|) must vanish 
if the degeneracy factor q > 2. This is indeed the case as 
for q > 1 = for at least one value of m ^ n. 



E. The DFT sum rules 

The HFB sum rules derived in Sec. Ill CI are based on 
the linearity of the Hamiltonian, by which a matrix ele- 
ment involving the HFB state is a sum of matrix elements 
calculated for all the PNP components (fl~8|) . In order to 
derive the analogous sum rules for the projected DFT 
energies, one can only use properties of the underlying 
transition energy density. To this end, we recall that 
in the HFB theory, the mixing of particle numbers corre- 
sponds to the broken U(l) gauge symmetry, and that the 
PNP actually corresponds to expanding the HFB state in 
irreducible representations of this group. This observa- 
tion can be extended to the DFT transition energy den- 
sity, expanded in these same irreducible representations, 
with the projected DFT energies being the expansion co- 
efficients. The resulting sum rules must follow from the 
closure relations on the group manifold. 

These general remarks can be expressed in an explicit 
form in the following way. By using integration contours 
that are circles of radius \zq\ around the origin, z = z^e % ^ , 
we have the following expression for the projected DFT 



and C p being a coupling constant. Assuming a two- 
fold Kramers degeneracy, the corresponding residue at 
±i\u n /v n \ is: 



(47) 

I 

energy 

-N f 2it 

(* N \* N )E» FT = / d<l>e- iN <t>£(ct>), (48) 

27r Jo 

where by £{<$>) we denoted the part of the integrand that 
does not depend on N, i.e., 

8(0) = (<S>Mz))E BFT (p z , Xz ,Xz) at z = z e i *. (49) 

Hence, the DFT projected energy is given by a Fourier 
transform of £(<fi). Since the Fourier components 
constitute a complete set of functions on a circle, 
J2n=o e ~ m,p = 2tt(5(0), we obtain the DFT sum rule, 

OO 

($\<S>(z ))E DFT ( Pzo , Xzo ,Xzo) = E Zo(*n\*n)E$ ft , 

which is the analogue of the HFB sum rule for matrix 
elements (|24| . For zq—1, we obtain the DFT counterpart 
of the HFB sum rule 

OO 

E DFT (p, X ,X*) = E^I^Xft- (51) 

N=0 

We note that in the above derivations, zq is an arbi- 
trary complex number; its modulus fixes the radius of 
integration contour, while its phase gives the point on 
the circle that fixes the starting point of the integral in 
Eq. (|48[) . This starting point has obviously no impor- 
tance for the value of the integral. The sum rule ([50]) 
gives, therefore, a representation of the DFT transition 
energy density in terms of a series expansion in zq, which 
converges only on the ring between the poles. For each 
such ring, the projected DFT energies -S'dpt are differ- 
ent, and the DFT transition energy density is thus equal 
to a different series expansion. It is obvious that these 
different values of the projected DFT energies do not 
contradict the continuity of the DFT transition energy 
density. In this way, all projected DFT energies for arbi- 
trarily chosen contours of integration correspond to this 
same common DFT energy functional. 



n a 

I 



, . - (-1) " I* (e i*H'.n < (f - 1) 



F. Density-dependent terms with fractional powers 

Let us now analyze the terms in the DFT energy den- 
sity that depend on fractional powers a of the local den- 
sity. In many functionals related to the Skyrme interac- 
tion, and for the Gogny force, such terms are quite often 
postulated, both in the particle-hole and pairing chan- 
nels (see Ref. [l2| for a review). In particular, the famil- 
iar density-dependent term of the Skyrme force, which is 
proportional to p 1 (r) , produces a contribution of the or- 
der of a = 2+7 to the DFT energy density. Similarly, the 
density-dependent, zero-range term of the Gogny force 
yields a contribution to the DFT energy density that is 
of a — 1 + 7 order, provided the particle-hole and pair- 
ing terms are consistently added. By taking into account 
the degeneracy factors q discussed in Sec. IIIICl the re- 
sulting poles are of the order of a — p = 2 + 7 — p and 
a — p = 1 + 7 — p, respectively. Since typical values of 7 
are between and 1, the DFT transition energy density 
always has poles at z n = ±i\u n /v n \ for non-degenerate 
states (q — 1). But more importantly, the fractional pow- 
ers lead to the multivalued DFT transition energy density 
on the complex- z plane. 

In the standard treatment of fractional powers a of a 
complex function, cuts along the negative real axis must 
be introduced. In order to apply this procedure to frac- 
tional powers of the local transition density (I46|) , we must 
identify on the complex z plane the lines along which 
p z (r) is real negative. 

Obviously, p z (r) is real positive along the real z axis 
and real along the imaginary z axis. In order to simplify 
the discussion, let us assume that the sum in Eq. (|46) is 
finite, which is always the case in any practical calcula- 
tion. In such a case, p z (v) has a finite number, say M, 
of different first-order poles along the positive imaginary 
axis, and the same number M of poles along the neg- 
ative imaginary axis. Moreover, since all coefficients in 
Eq. (|4"6"|) are positive, p z (r) must have a first-order zero 
between each pair of poles on the positive imaginary axis, 
and similarly on the negative imaginary axis. Since p z (r) 
also has a second-order zero at z = 0, we conclude that it 
has altogether 2M zeros on the imaginary axis. It is also 
obvious that p z (r) is a rational function with a 2M-order 
polynomial in the numerator, and thus we conclude that 
all the zeros of p z (r) are located on the imaginary axis. 
Therefore, the cuts for possible fractional powers a must 
be located along the imaginary axis, and connect zeros 
of p z (r) with its adjacent poles. 

The above discussion is visualized in Fig. O The 
left portion shows schematically the transition density 
Pim[z]( r ) along the imaginary axis Im[z] oriented verti- 
cally. The plot illustrates the transition density (|4"6"| in 
one selected point of space r, i.e., values of wave functions 
at r enter only as numerical coefficients. There appear 
four poles and three zeros of pi m \ z ] on the positive imag- 
inary axis, the same number of poles and zeros on the 
negative imaginary axis, and the second-order zero at the 
origin. Sections of the imaginary axis where the density 
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FIG. 3: (color online) Ordinary (2 = 1) and transition (z = i) 
densities in 18 calculated in HFB+SLy4 as functions of r. 
The upper part of the Figure shows (in extended scale) the 
small region of radii where the transition density becomes 
negative. 

is negative are shaded. In the right portion of Fig. [2] we 
show poles (crossed circles) and zeros (full dots) of the 
transition density on the complex z plane, along with the 
three integration contours CO, CI', and C2' discussed 
above. The cuts in the complex z-plane connecting zeros 
and poles, corresponding to real negative values of p z , 
are indicated by vertical segments. 

While the location of the poles is independent of r, the 
position of zeros of the transition density is r-dependent. 
In order to visualize this, in Fig. [3] we plot the total or- 
dinary (z — I) and transition (z = i) densities in 18 
obtained within the SLy4 energy density functional. One 
can see that the transition density is positive almost ev- 
erywhere; only in a very narrow region near r ~ 5 fm 
does it become slightly negative, as shown in the upper 
part where the scale is expanded by a factor of 1000. 
Such a behavior of p z {r) at z 2 = —1 can be easily un- 
derstood from Eq. (|46|) . Indeed, strongly occupied states 
with ~ I always yield positive contributions, while 
negative contributions of states with w~ < can only 
appear in the surface region where the least bound canon- 
ical states dominate. 

By the same token, we can see that strong negative 
contributions may appear when the integration contour 
passes slightly below a pole located just above the unit 
radius, i.e., is slightly smaller then u\. Such a situa- 
tion is predicted in 26 0, where the occupation probability 
of the canonical 2d 3 / 2 state equals 0.486. As shown in 
the bottom part of Fig. [4j the transition density at z = i 
is dominated by this particular contribution and becomes 
strongly negative beyond r ~ 2 fm. 

We are now ready to discuss contour integration of 
terms depending on fractional powers of the transition 
density. Contour CO shown in Fig. [2] crosses the imagi- 
nary axis in sections where there is no cut, and thus it 
always stays on the same Riemann sheet. On the other 
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FIG. 4: (color online) Bottom: ordinary (z = 1) and tran- 
sition (z = i) densities in 2e O calculated in HFB+SLy4 as 
functions of r. Top: the real parts (solid lines) and imaginary 
parts (dashed lines) of the 1/6 powers of the ordinary (z = 1) 
and transition (z = ie ±ts ) densities. 



hand, contours CI and C2 of Fig. [TJ or contours CI' and 
C2' of Fig.[2l cross the imaginary axis by passing through 
cuts onto another Riemann sheet. Since the transition 
density ([4^)1 is an even function of z, the phase of the 
fractional power a of the transition density increases or 
decreases by 2ira when going across each of the two cuts. 
Therefore, after returning to z = 1, p z (r) is multiplied 
by exp(±47ra), and thus it is not a continuous function 
at z = 1, unless a = k/2. This is quite unacceptable 
as the presence of the phase creates serious problems in 
interpreting the projected DFT energies (see, e.g., 
the sum rule condition discussed in Sec. IIII E|) . 

Formally, by using powers of square roots in density- 
dependent terms, i.e., a = k/2, one can guarantee that 
the integration contours return onto the original Rie- 
mann sheet, and that the transition energy density is 
a continuous function of z. However, even in such a case, 
one important property of the DFT transition energy 
density (|33|) is lost, namely, the density-dependent term 
in the energy density becomes an odd function of z, and 





FIG. 5: Modified closed contour CI" that encircles the zero 
of the transition density. The poles (zeros) of p z are marked 
by crossed circles) (dots). The equivalent circular contour 
CI'" lying between the zero of p z and the previous pole is 
also indicated. See text for more details. 



the corresponding term in the integrand (|35j) becomes an 
even function of z. This is so, because the square root 
has opposite signs on the two Riemann sheets in ques- 
tion. Consequently, contour integrals of such terms would 
vanish and the density-dependent terms would yield zero 
contribution to the PN-projected energy. This is a rather 
disastrous result. Hence, we are forced to conclude that 
the use of continuous contours for fractional powers is 
not a viable prescription for constructing the projected 
DFT energies. 

Let us now discuss the way of evaluating contour inte- 
grals in all practical PNP calculations up to now. Unfor- 
tunately, such calculations have always disregarded the 
analytic structure of the underlying integrands. In fact, 
the fractional powers of transition density, 



P™( r ) = M r )|exp{iarg[p z (r)]/a}, 



(52) 



are practically determined by computer compilers. In 
Eq. (p)2"j). the so-called argument arg [p z (r)] of p z (r) is de- 
fined as the phase of the complex variable p z (r); hence, it 
is contained in the interval of — tt^tTt. This usual prescrip- 
tion corresponds to stepping over the cut whenever the 
contour approaches the imaginary axis for arg [p z (r)] = 
±7r, i.e., for real negative transition densities Pz(r) < 0. 
In this way, the integrand is always calculated on the 
same Riemann sheet, but the integration contour is not 
closed. 

The contour can be closed by adding a piece that goes 
around the zero of the transition density p 2 (r). This is 
illustrated in Fig. [5] that shows a modification of con- 
tour CI' of Fig. [5] near the positive imaginary z-axis. 
[An analogous mirror-like detour is made near the nega- 
tive imaginary axis.] The resulting contour CI" always 
stays on the same Riemann sheet and, therefore, the in- 
tegration result does not depend on the radius. On the 
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other hand, the contribution due to the additional path 
surrounding the zero is affected by the discontinuity of 
the integrand along the cut. Such a discontinuity in the 
transition density is shown in the top panel of Fig. 2] for 
a=l/6 ([52]) . Since the ordinary density (z — 1) is real 
and positive, its fractional power is also real and positive. 
On the other hand, transition densities (|46|) correspond- 
ing to z± = ie ±i£ with e — 0.5°, i.e., near the positive 
imaginary axis, are complex. While their real parts are 
practically identical on both sides of the cut, their imag- 
inary parts have opposite signs; hence, a discontinuity is 
encountered. [Of course, since the directions of integra- 
tion are opposite, the contributions to the PNP energy 
from both segments z± of the additional path are identi- 
cal.] 

When transforming the contour integration in Eq. (1341) 
into the integral along the imaginary axis y, one must do 
a change of variables from z to iy. This introduces the 
additional factor i~ N in the integrand. For that reason, 
for even N, the discontinuity in the imaginary part of the 
density-dependent term of fractional order contributes 
to the real part of the projected DFT energy. As the 
discussion in Sec. IIII Dl proves, the same holds for odd 
values of N. 

Figure[5]also shows the circular contour CI"' lying be- 
tween the zero of p z and the previous pole |z„_i|. This 
contour is formally equivalent to the deformed contour 
CI" but it is easier to handle in practical applications. 
The radius of CI"' must be slightly greater than |z„_i| 
and smaller than the lowest zero of p z (r), minimized over 
the whole space r, associated with the branching point 
corresponding to z n . The use of contour CI'" guaran- 
tees that the integration of fractional-order terms is done 
properly. 

Altogether, blind application of prescription ([52")) can 
lead to spurious and entirely uncontrolled contributions 
to the projected DFT energies. Excepting Ref. [23j, this 
fact has been entirely overlooked in all practical appli- 
cations of the PNP method to date, and casts serious 
doubts on the reliability of the obtained results. Largest 
contributions are, of course, obtained when the integra- 
tion contour passes slightly below a pole of the DFT 
transition energy density. For the Skyrme functionals in 
Sec. UVBl we present specific examples of such situations. 

The appearance of spurious contributions is, in fact, in- 
dependent of the order of divergence at the pole. There- 
fore, it also shows up for "integrable" poles, diverging 
with powers of a — p = 1 + 7 — p < 1, discussed for the 
Gogny force in Ref. [37| . 



IV. NUMERICAL EXAMPLES 

In order to illustrate theoretical findings presented in 
Sec. IIII1 we carried out numerical calculations within the 
Skyrme-DFT method. We used the code HFBTHO [H 
which is capable of handling spherical and axially de- 
formed nuclei within the Lipkin-Nogami (LN) approx- 



imation followed by the PNP. This corresponds to the 
projection-after-variation (PAV) method of restoring the 
PN symmetry. By using a new version of HFBTHO, we 
also performed full variation-after-projection (VAP) cal- 
culations analogous to those of Ref. [23|] . 

As illustrative examples, we study spherical and de- 
formed configurations in 18 O and in 32 Mg calculated us- 
ing the Skyrme functionals SIII [Io[ and SLy4 [|l| . These 
two parametrizations differ in a significant way with re- 
spect to the PNP method. The density-dependent term 
of SIII contributes to the energy density as (p n +Pp)pnPp- 
Therefore, both in the neutron and proton subsystems, 
the powers p (Sec. IIII Ep of the density dependence are 
equal to 2. Consequently, from the PNP perspective, 
the density-dependent term of SIII is not any different 
than the density-independent terms. On the contrary, 
the density-dependent term of SLy4 is proportional to 
[Pn + Pp] 1 ^ 6 and exemplifies the case of fractional-power 
dependence discussed in Sec. IIII Fl The contact pairing 
force of the volume type (density-independent) was used 
in the particle-particle channel. All calculations have 
been performed in the spherical harmonic-oscillator basis 
of Nq = 6 or 10 shells, for 18 or 32 Mg, respectively. 



A. Numerical accuracy 

To calculate residues, we take circular contour integrals 
of radius ro : 

z = r e^. (53) 

The integrals are evaluated using the Fomenko discretiza- 
tion method [42|, E^|, whereby values of integrands are 
summed up at gauge angles <j>k = %r for k = 0, . . . , L— 1. 
This corresponds to the upper half circle in the com- 
plex z-plane and, as discussed in Sec. IIII Bl only the real 
part of the integral is kept. For analytic integrands, the 
Fomenko method delivers exact results up to admixtures 
of wave functions with N ± L, N ± 2L, N ± 3L, . . . par- 
ticles. The main question in applying this method to 
non-analytic integrands, which have poles in the complex 
plane, is to what extend can it deliver equally accurate 
numerical results. 

The Fomenko method clearly fails when there is a pole 
(|28[) lying just on the integration contour, ro = \z n \, and 
an even number of points L is used. In such a case, 
the integration point with k=L/2 is located exactly at 
the pole of the integrand. Therefore, in most practical 
calculations, an odd number of integration points, most 
often L = 7 or 9, was used. 

However, a more stringent condition on L results from 
the fact that the discretization method must fail when- 
ever the integrand varies too rapidly between two neigh- 
boring integration points. Therefore, the spacing be- 
tween points irro/L must be appropriately smaller than 
the distance from the pole. For odd values of L, the inte- 
gration points corresponding to k = (L± l)/2 are closest 
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to the imaginary axis; hence, one arrives at the condition 



7rr 



L <Vl"»-W>* + (^). 



L > 



V^\r - \z v 



(54) 



(55) 



In the present study, a large number of L — 93 integration 
points was used, which allows for calculating the contour 
integrals with radii ro that differ by as little as 3% from 
the position of the closest pole \z n \. 



B. Dependence of projected energy on integration 
contours in spherical nuclei 

Table fl] displays the results of PNP calculations per- 
formed for 18 O by using circular integration contours (|53[) 
of different radii. The precision of numerical integrations 
was confirmed by calculating contributions from individ- 
ual poles. This was done by carrying out contour inte- 
grals over small circles surrounding the poles. In this 
way, we determined residues from the individual poles 
E^ FT (n) and checked that their sums, X^m=o -^dft ( m ) > 
agree very well with the results of contour integrals along 
circular contours C, as required by the Cauchy theorem 

da. 

As seen in Table HI contributions of the n = poles at 
z = are huge. Therefore, the DFT residues at z = 
cannot at all be interpreted as the projected energies, as 
was the case for the PNP HFB theory, Eq. (|15p. Residues 
at z — are cancelled, to a large extent, by contribu- 
tions from the lsi/2 deep-hole states, which are large be- 
cause they contain large factors of the type (— v^/u^) 
for u„ ~ [see Eq. (j47] . Contributions from other poles 
are also quite large, and apart from the integration con- 
tour at \z\ = 1, none of the other contours reproduce the 
correct projected energy shown by a boxed number. 

For the SIII parametrization, one can see that contri- 
butions from poles associated with spherical states with 
j > 3/2 (g > 2) are indeed equal to zero, cf. discussion 
in Sec. IIII CI This property does not hold for SLy4, for 
which the projected DFT energies have jumps also when 
the integration contours cross the j > 3/2 poles. In this 
case, the jumps are not related to non-zero residues, but, 
as discussed in Sec. IIII Fl they are caused by the fact that 
the integration contours are not closed for the fractional- 
power terms. 

Figure [5] shows the projected DFT-SLy4 energies ob- 
tained by using circular integration contours of different 
radii ro. These calculations illustrate properties of poles 
listed in Table |TJ The contributions originating from the 
density-independent and density-dependent terms of the 
Skyrme force are separated. The latter terms yield the 
fractional-power terms in the DFT energy density dis- 
cussed in Sec. IIII Fl As in the SIII case, the density- 
independent terms exhibit jumps only at the two j = 1/2 
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FIG. 6: (color online) The N = 10 projected DFT-SLy4 en- 
ergies (|44p calculated in 18 O as functions of the integration 
radius ro. Upper (a) and lower (b) panels show, respectively, 
results for the terms originating from density-independent and 
density-dependent parts of the Skyrme force. Positions of in- 
dividual neutron poles are marked by arrows. The inset shows 
the results near the 2si/2 pole on an expanded scale. The 
result of calculations obtained with the equivalent contours 
passing below the branching point associated with the 2si/2 
pole (cf. Fig. [5} is shown by a dotted line. 



poles. On the other hand, the density-dependent terms 
show jumps at all poles, and these jumps carry over to 
the total projected DFT energies shown in Table fl] [The 
small jump at the lds/2 pole, 120 keV, is practically in- 
visible in the scale of Fig. Moreover, contributions 
of the density-dependent terms are not constant between 
the poles, as would be required by the Cauchy theorem. 
This is caused by the prescription (|52|) to step over the 
cuts in the complex plane, and illustrates spurious con- 
tributions to the projected DFT energies discussed in 
Sec. IIII Fl As shown in the blown-up inset in Fig. [5Jb), 
these spurious contributions appear just below the pole 
thresholds (i.e., for small negative values of ro — \z n \), 
and they can be quite large - of the order of several tens 
of MeV. The gradual development of spurious contribu- 
tions below threshold has been explained in Sec. IIII Fl 
Namely, if the contour radius is only slightly greater than 
|z„_i|, the branching point associated with the pole z n 
is always outside for all values of r. With increasing ro, 
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TABLE I: Contributions E^ FT (n) in 18 O from the individual neutron poles to the projected DFT energy (|44|) for N — 10 
calculated using the SIII and SLy4 Skyrme functionals. For each parametrization, the last column shows the sum of the nth 
lowest contributions (144 \ , with the values of S DFT marked by boxed numbers. Canonical energies e n and pole positions z n are 
also given. All energies are in MeV. 





SIII 


SLy4 


n orbital 




£n Zn -EdFtH Em=0 "^DFT (w) 


n.a. 

1 lSl/2 

2 IP3/2 

3 lPl/2 

4 ld 5/2 

5 2si/ 2 

6 ld 3/2 

7 lf 7 /2 


n.a. 0.000 -1.5io+6 -1.5i +6 
-35.181 0.021 1.5io+6 0.178 
-20.302 0.045 0.178 


n.a. 0.000 -2.9io+6 -2.9i +6 
-36.985 0.038 2.9 M +6 731.008 
-20.691 0.077 -8.6io+2 -125.448 


-15.040 0.072 -1.4k,+2 -140.540 


-14.783 0.131 -1.7io+2 -142.584 


-6.528 1.429 -140.540 
-2.166 11.255 -l.lio+3 -1235.194 
2.831 17.458 -1235.194 
10.265 31.398 -1235.194 


-6.399 1.462 1.2i -l -142.464 
-2.685 5.702 -4.8 UJ +1 -190.628 
3.143 10.952 4.0io+l -150.627 
10.325 19.033 2.5io+l -125.608 



more and more branching points corresponding to differ- 
ent regions of space fall inside the contour, leading to the 
spurious behavior. As discussed earlier, one can elim- 
inate this subthreshold effect by taking equivalent con- 
tours discussed in the context of Fig. [5j Such a procedure 
is illustrated by a dotted line in the inset of Fig. [6fb). 

The spurious contributions may result in large errors 
in the projected PNP energies, making results of the 
standard PNP calculations meaningless. Unfortunately, 
this is true not only for Skyrme forces that use density- 
dependent terms of fractional orders but also for the 
Gogny force, which contains a density-dependent term 
of order 7 = 1/3. 



C. Calculations for deformed nuclei 

In our previous study [44J, we calculated the com- 
plete HFB mass chart of even-even nuclei by performing 
the PNP of paired ground states determined by the LN 
method. At this point, when performing the PNP calcu- 
lation in each individual nucleus, one should take care of 
the cases when one of the poles z n turns out to be near 
the standard tq — 1 integration circle (unit circle). 

In order to produce the ground state masses for all 
even-even nuclei lying between the two-nucleon drip lines, 
one has to calculate about 6000 nuclei. Moreover, each 
nucleus has to be calculated three times, by starting from 
oblate, spherical, and prolate initial shapes. We have 
found that among these 6000 nuclei, about 100 have a 
neutron or proton state with occupation numbers near 
1/2. Therefore, the standard PNP method yields about 
100 questionable results across the mass chart. However, 
the situations is much more serious when performing the 
constrained HFB calculations discussed in the following 
sections. 



D. Distribution of poles as a function of 
deformation 

When increasing the quadruple deformation, states 
with smallest (largest) angular momentum projections 
onto the symmetry axis, f2, become more (less) bound 
on the prolate side, and the opposite holds for the oblate 
side. For states located above the shell gap, this means 
that low-fl and high-Si orbitals become more occupied 
with increasing prolate and oblate deformation, respec- 
tively. Therefore, at some deformation, these orbitals 
cross the Fermi energy and the corresponding poles cross 
the unit circle. An analogous situation may also occur 
for orbitals located below the shell gap, whereupon high- 

and low-fi Nilsson orbitals become less occupied with 
increasing prolate and oblate deformation, respectively, 
and also may cross the Fermi energy. We wish to empha- 
size that the problem occurs not at the point where the 
orbitals from above and below the shell gap cross each 
other, leading to a configuration change, but at defor- 
mation where either of these orbitals crosses the Fermi 
energy. 

Such a case is illustrated in Fig. [7] for the nucleus 18 0. 
In pure HFB calculations (no LN correlations included) , 
this nucleus has neutron pairing only. At the spherical 
shape, the ld 5 / 2 shell is located above the A^=8 shell 
gap, i.e., it has particle character (|z| > 1). The three- 
fold degeneracy of this shell (g=3) makes the contribu- 
tion from this pole to the projected energy vanish. At 
nonzero deformations, however, the degeneracy is lifted 
and three individual poles (<?=1) appear in the complex 
plane. Moreover, near (3 = 0.12 and (3 = —0.12, poles 
corresponding to the 0=1 /2 and $1=5/2 Nilsson levels 
cross the unit circle \z\ = 1. 

The situation is much worse for nuclei having more 
single-particle states with poles close to the unit circle. 
The neutron- rich nucleus 32 Mg is such a complicated case 
illustrated in Fig. [5] This example is calculated in the 
HFB+LN approach, in which both neutron and proton 
pairing is nonzero. For completeness, canonical single- 



13 



100 
10 

c 

N 1 

0.1 
0.01 
0.001 



wvvvvvvvvvv vv w w 

2s. „ 



1d 



.nnn n[][]D D D □ n n 

^□□dJuff^ 

1d__ 1d.„.„ ■■■ 



1P 

1P, 



^^^YaWWAAAAAAA A A A A A A A A AAAAAAAAAAAA 



O * O ^ O $ & * 



"IS AAAAAAAWAAAAAAAAAAAAAAAAAAAAAAA a a , 

1/2 ""AAAAAaa 



-0.2 -0.1 0.0 0.1 0.2 
deformation p 

FIG. 7: (color online) Neutron poles z„ in 18 (|28[) as func- 
tions of quadrupole deformation /3, calculated within the 
HFB-SIII method with volume pairing interaction. 



particle energies e n associated with the poles z n are plot- 
ted in Fig. M 

As can be seen in Fig. [5J there appear numerous cross- 
ings of poles with the unit circle as a function of defor- 
mation. On the prolate side, neutron poles lf7/ 2 ,n=i/2 
(Nilsson level [330] 1/2) and ld 3A3 / 2 ([202]3/2) cross the 
unit circle at the same deformation where they cross one 
another. At larger deformation, the same situation oc- 
curs for the lf 7/2 , 3 / 2 Q321]3/2) and ld 3/2 , 1/2 ([200]l/2) 
orbitals. For protons, a single ld 5 / 2 5/ 2 ([202] 5/2) or- 
bital crosses the unit circle at small deformations. On 
the oblate side, neutron orbitals lfV/2,5/2 ([312]5/2) and 
ld.3/2,1/2 cross the unit circle at different deformations, 
near the point where they cross one another, while the 
proton ld 3 / 2 ,i/2 orbital stays near the unit circle for a 
wide range of deformations. 

As discussed in Sees. IIII CI and IIVB1 results of the 
PNP, at least for the density- independent terms, must 
only depend on the residues of poles that are inside 
the integration radius r$. However, whenever a given 
pole crosses the integration contour, the projected en- 
ergy must undergo a sudden jump as a function of defor- 
mation. This jump is, of course, equal to the residue at 
this pole. The fact that a given pole crosses the integra- 
tion contour could be without consequence, provided the 
contour is shifted back to always stay between the same 
poles. This is always possible, as long as the poles do 
not cross around the contour. It is obvious that when- 
ever they do, the projected energy may have a sudden 
jump that cannot be avoided by a contour shift. On the 
other hand, when two poles cross precisely at the inte- 
gration contour, the corresponding degeneracy factor q 
increases by a unity, and the poles may simply disappear 
(at least for the terms that show polynomial density de- 
pendence), in which case the projected energy may stay 
smooth. Such cases are studied in the next section. 



E. Deformation energy within the HFB+PNP 
method 



Results in this section were obtained with the SIII 
Skyrme force whose density-dependent terms do not cre- 
ate additional problems (cf. Sec. IIV[) . Let us first ana- 
lyze the simpler case of 18 O. Figure 1X01 presents the de- 
formation energy E((3) as a function of the quadrupole 
deformation (3. As a reference curve, we show the un- 
projected deformation energy emerging from the HFB 
calculations and the associated PNP energy curve (PAV; 
solid squares; the contour radius ro=l). Near the ground 
state ((3 — 0), the projected energy is lowered by about 
1.5 MeV due to additional PN correlations. At larger de- 
formations, the correlation energy decreases due to the 
stronger static neutron pairing. 

At deformations ~ —0.12 and (3 ~ +0.12, the pro- 
jected energy curve exhibits unphysical jumps. By com- 
paring with Fig. one concludes that at these defor- 
mations the neutron ld5/ 2 poles cross the integration 
contour. Obviously, the residue contributions of these 
poles cause the sudden jumps in the deformation energy. 
The ld 5 / 2j 5/ 2 pole introduces a positive contribution at 
(3 ~ —0.12, while the ld5/ 2j i/ 2 pole introduces another 
positive contribution at (3 ~ +0.12. Based on this obser- 
vation, two other sets of PAV calculations were carried 
out. The first calculation (open circles) was done by ex- 
cluding contributions from the ld5/ 2 poles, as is the case 
for the ground state configuration. This can be accom- 
plished by reducing the integration radius from ro = 1 
to a smaller value of about ro = 0.1, cf. Fig.[7J At small 
deformations, —0.12 < (3 < 0.12, the new results are 
identical to those obtained with the unit circle, and at 
the larger deformations (prolate or oblate), the energy 
curve smoothly continues without any jump. Thus, in 
this example, an appropriate shift of the integration con- 
tour allows us to obtain smooth and unique projected 
energy. The second PAV curve (open squares) has been 
obtained by always including the lowest lds/ 2 pole, i.e., 
by continuously varying the integration radius as a func- 
tion of (3, to ensure that it always stays between the first 
and second ld 5 / 2 pole (cf. Fig. UJ. The resulting energy 
curve coincides at large deformations with the standard 
PAV result, and then smoothly continues to (3 = 0, where 
the ld 5 / 2 poles (q=3) disappear. 

Figure [TO] also presents the fully self-consistent VAP 
results. Similar to the PAV case, two sets of calculations 
were performed. The solid (open) stars correspond to 
including (excluding) the contributions from the lowest 
Ids/2 poles. In both cases, one obtains smooth curves, 
which, beyond the spherical point, differ from one an- 
other. 

In this rather simple case of 18 O, both in the PAV 
and VAP calculations, one can avoid unphysical jumps 
of the projected energy curve by making a specific selec- 
tion of "active" poles that are considered during contour 
integration. Such selection of residues can, in principle, 
become a part of the definition of the projected energy. 
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FIG. 8: (color online) Neutron (left) and proton (right) poles z n (|28l) as functions of quadrupole deformation /3 calculated for 
32 Mg with the SIII functional and volume pairing interaction. 
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FIG. 9: (color online) Similar to Fig. [8] except for canonical energies e^. 



The variational principle can then be invoked to pick the 
selection that yields the lowest projected energy. In the 
discussed case of 18 O, the PAV and VAP energies ob- 
tained by excluding the lds/2 poles are the lowest, and 
they are smooth functions of deformation. Therefore, 
such a selection can be adopted for the final PNP energy 
in this nucleus. It is clear, however, that one cannot a 
priori tell which selection of poles leads to the lowest pro- 
jected energy. For example, in heavier oxygen isotopes, 
the lowest energy is obtained by including some of the 
ld 5 / 2 poles. 



Let us now consider a more complicated case of the 
HFB+LN calculations for the neutron rich nucleus 32 Mg. 
The total HFB energy (without the corrective A^ 2 -* LN 
term) is shown in Fig. [11] as a function of f3. Solid squares 
denote the result of PAV PNP calculations on the top of 
HFB+LN. At (3 ~ +0.1, the PAV curve exhibits a small 
jump, after which its behavior changes character. This is 
clearly related to the proton ^5/2,5/2 pole crossing the 
unit circle, cf. Fig. [51 Otherwise, the PAV deformation 
energy is quite smooth as a function of deformation, de- 
spite the fact that three pairs of neutron poles cross the 
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FIG. 10: (color online) Deformation energy E(f3) as a func- 
tion of quadrupole deformation /3, calculated for ls O with the 
SIII force and volume pairing interaction. Results of HFB 
(open triangles) are compared with different variants of PAV 
(squares and circles) and VAP PNP (stars) (see text for de- 
tails). 



unit circle in the deformation range considered. This ap- 
parent lack of sensitivity to neutron poles can be traced 
back to the fact that they cross the integration contour 
at or near points where they pairwise cross one another. 
Therefore, the increasing degeneracy factor q makes the 
poles disappear at the crossing points; hence, the total 
PAV curve behaves smoothly. 
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FIG. 11: (color online) Deformation energy E(f3) as a function 
of quadrupole deformation f3 calculated for 32 Mg with the 
SIII force and volume pairing interaction. Results of the PAV 
HFB+LN calculations (squares and triangles) are compared 
with the VAP PNP results (dots). The standard HFB result 
is shown by open triangles. 

Following the example of 18 O, also for 32 Mg we per- 
formed PAV calculations wherein we took into account 



contributions from all the poles that contribute to the 
ground state configuration (/? = 0). The resulting PAV 
curve (solid triangles in Fig. ITT]) behaves smoothly but 
shows unphysical behavior at very large deformations. 
An explanation of this artifact follows from Fig. [T2l 
where we plot energy contributions from the most im- 
portant poles: (i) the neutron ld 3 / 2 ,3/2 pole (solid cir- 
cles; it leaves the unit-circle at /3 ~ 0.22 and we have 
to add its contribution beyond this point); (ii) the neu- 
tron ld7/2,i/2 pole (solid squares; it enters the contour 
at /3 ~ 0.22 and we have to subtract its contribution be- 
yond this point); and (hi) the proton ^5/2,5/2 pole (solid 
triangles; it leaves the contour at (3 ~ 0.1 and we have to 
add its contribution beyond this point). Interestingly, 
pole contributions to the total projected energy oscil- 
late with deformation. As expected, the residues van- 
ish when the corresponding poles cross at the integration 
circle, cf. Fig.[H] However, oscillations between the cross- 
ing points can become quite large, as is the case for the 
proton lds/2,5/2 pole; hence, the projected energy curve 
obtained by keeping contributions from the ground-state 
set of poles acquires strong unphysical oscillations. 
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FIG. 12: (color online) Contributions to the PAV energy of 
32 Mg from the three selected poles as a function of deforma- 
tion 13. 

In a search for the most sensible method of calculating 
the projected energy curve within the PAV approach, we 
can employ a prescription whereby the number of the 
lowest poles is kept fixed within the integration contour. 
This can be realized by keeping vq between the poles or at 
the pole crossing point. Since at the crossing points the 
poles vanish, at least in SIII calculations, this results in a 
smooth energy curve. Such an option is shown in Fig. 1111 
with solid squares. The resulting curve is indeed very 
smooth; however, at /3 ~ 0.4 there appears an unphysical 
bump, which makes this option as unacceptable as the 
other one. 

We have also attempted calculating the energy curve 



16 



within the VAP approach. In principle, the VAP ap- 
proach could have generated problems related to the fact 
that the density-dependent terms of fractional order may 
lead to large negative contributions (see Sec. IIII F|) . In 
practice, this is never the case because the VAP method 
[23j is not implemented through an explicit minimization 
of the energy, but is carried out by solving variational 
equations that have been derived with the same incorrect 
treatment of cuts in the complex plane. In this context, it 
is worth emphasizing that the appearance of poles never 
leads to infinite total energies, but to discontinuities in 
the total energy. Therefore, there is no danger that the 
minimization procedure may attract a solution towards 
a pole. 

The main problem in implementing the VAP method 
is related to the fact that unprojected quantities, e.g., 
particle p nn > or pairing p nn i densities, lose their usual 
physical meaning [23j in VAP. They depend on the in- 
ternal normalization N' = Trp of the density p nn / that 
is not related to the particle number N onto which the 
state is projected. Neither the total VAP energy nor 
other projected observables depend on the normalization 
N'. However, depending on the choice of the internal 
normalization N', one obtains different canonical occu- 
pation probabilities; hence, the associated poles z n are 
not distributed in the same way as in the unprojected 
HFB case. Depending on the internal normalization N' , 
different poles z n enter the integration contour, and the 
convergence procedure cannot be easily controlled. 

Additional problems arise when two poles are nearly 
degenerate. Although at the point of degeneracy the 
poles disappear, when the distance between the poles is 
small but nonzero and the integration contour is between 
them, one faces significant instabilities of the constrained 
VAP problem. During the iteration of VAP equations, 
one or both poles enter or leave the integration con- 
tour. The poles create jumps in the projected energy 
and, which is even more important, they create jumps in 
the deformation. The numerical algorithm enters a 'ping- 
pong' regime, which cannot be overcome, and one cannot 
converge to any solution. Figure [IT] offers a good illus- 
tration of this problem. The converged VAP energies for 
32 Mg are shown with solid circles. The converged solu- 
tion can be found only in limited regions of deformation. 
The Id and If neutron poles close to the contour spoil 
the convergence in the regions of (3 ss —0.2, 0.25, and 
0.5. The same is true for the proton Id states around 
j3 ~ —0.2 and 0.1. As a result, the VAP procedure could 
be solved only within small deformation intervals around 
(3 w -0.2, 0, and 0.4. 



V. CONCLUSIONS 

This study contains the systematic analysis of the 
particle-number-projected DFT approach. This ap- 
proach, usually in the Skyrme-HFB or Gogny-HFB 
framework, is commonly used in systematic calculations 



of nuclear ground-state properties, low-energy excita- 
tions, and high-spin states. For heavy, complex nuclei, 
the nuclear DFT is the only viable microscopic tool based 
on the effective interaction (or functional). To advance 
the nuclear DFT further, to improve its theoretical foun- 
dations, and to make it a more reliable tool, it is impor- 
tant to fully understand the advantages and drawbacks of 
the method when applied to self-bound nuclear systems. 

The main conclusions of this work can be summarized 
as follows: 

1. The transition density matrices connecting states 
having different orientation in the gauge plane z 
have poles on the imaginary axis Im[z]. In the HFB 
formalism that is based on one Hamiltonian act- 
ing in all channels (Hamiltonian-based HFB), these 
poles are irrelevant as their impact is nullified by 
the cancellation between the Hamiltonian matrix 
elements originating from particle-hole and pairing 
channels. Such a cancellation is not present in the 
HFB applications in which some of the Hamiltonian 
matrix elements are neglected (or approximated), 
in the HFB method based on density-dependent 
interactions (usually acting in different channels), 
and in the DFT approach in which the Hamilto- 
nian does not appear at all. In all those cases, the 
projection operator is not defined uniquely and the 
result depends on the analytic structure of the tran- 
sition energy density; hence, the projected DFT en- 
ergy depends on the integration contour. The re- 
sulting PNP energy can be expressed in terms of in- 
dividual residues corresponding to the poles z n as- 
sociated with single-particle (canonical) proton and 
neutron orbits. In contrast, in the Hamiltonian- 
based HFB, the result depends only on the single 
pole at the origin (zq — 0). 

2. Within the Hamiltonian-based HFB, there exist 
sum rules that relate the unprojected matrix ele- 
ments with the matrix elements in the projected 
states. Similar sum rules can be derived within the 
DFT; they relate the unprojected DFT transition 
energy density with the projected DFT energies. 
The DFT sum rules offer the interpretation of the 
projected DFT energies as Fourier components (as- 
sociated with irreducible representations of gauge 
group U(l)) of the DFT transition energy density. 
This can be naturally extended to other (higher) 
broken symmetry groups, such as SU(2) (associated 
with the broken angular momentum symmetry). 

3. The discussion of the particle number restoration 
can be extended to other symmetry restoration 
problems. In particular, DFT transition densi- 
ties associated with angular-momentum-projected 
states are expected to have complicated pole struc- 
ture in the three dimensional space of Euler angles 
(see the example shown in Ref. [45[). 

4. For the terms in the density functional that have 



17 



polynomial density dependence, the appearance of 
poles inside the contour gives rise to sudden jumps 
of the projected energy whenever the contour's pole 
content changes. Otherwise, the results are stable. 
This is not true for the terms having fractional- 
power density dependence (e.g., density-dependent 
pieces of many Skyrme and Gogny interactions or 
the Coulomb exchange term taken in the Slater ap- 
proximation). Here, the dependence on the contour 
radius shows a strong subthreshold behavior that 
can only be cured by considering appropriate inte- 
gration contours which do not go across the cuts in 
the complex z-plane. Other prescriptions give rise 
to uncontrolled energy behavior resulting from the 
fact that the corresponding integration contours do 
not close. 

5. As a practical measure that allows avoiding prob- 
lems related to the fractional-power density de- 
pendence, we propose using the integration con- 
tours which pass near and above the poles. Al- 
though such a prescription requires using rather 
dense meshes of integration points, it minimizes 
the risks of crossing the cuts in the complex plane. 
In this way, the ambiguities related to the non- 
analyticity of the DFT transition energy are re- 
duced to those corresponding to the choice of poles 
included within the integration contour. 

6. Projected DFT yields questionable results if a 
pole appears very close to the integration con- 
tour. While such a situation seldom happens in 
the ground-state calculations (less than 2% cases 
are affected) , it frequently occurs in calculations of 
projected energy surfaces, such as those in the gen- 
erator coordinate method (GCM). The appearance 
of poles in the vicinity of the contour as a func- 
tion of the collective coordinate (e.g., deformation) 
gives rise to uncontrolled irregularities and jumps 
in the results; in particular, it makes it impossible 
to define the PNP potential energy surfaces. 

7. Pole pathologies appear in a particularly strong 
way in the fully self-consistent VAP calculations. 
In this approach, transition density poles are not 
uniquely defined; moreover, their positions can 



change during the iteration process leading to nu- 
merical instabilities. 

8. The analytic structure of the transition energy den- 
sity becomes exceedingly complicated in nuclei with 
protons and neutrons paired, thus requiring simul- 
taneous proton and neutron PNP. Of particular im- 
portance in the context of GCM applications is the 
extension of the present analysis to non-diagonal 
matrix elements between the PNP states. 

Some of the problems listed above, in particular those 
related to the configuration-mixing DFT method and ap- 
plications of the generalized Wick's theorem to DFT, 
have been recently addressed in a series of papers [33l |34] | . 
In these studies, a practical cure has been proposed that 
is based on removing specific spurious components of 
the DFT functional that can be associated with self- 
interaction and self-pairing. When applied to truncated 
Hamiltonians, this practical prescription turns out to be 
very effective [4y]. However, this kind of solution does not 
remove ambiguities related to using complicated (e.g., 
fractional-power) dependence of the energy density func- 
tionals on particle densities. Finding ultimate cures to 
the problems discussed in our study will undoubtedly re- 
sult in establishing better theoretical constraints on the 
form of the DFT energy density functionals for nuclear 
self-bound systems. 
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